package fem2;

import fem2.enu.MeshPartType;
import fem2.mesh_and_topology.LineSpatialCriteria;
import fem2.mesh_and_topology.SpatialCriteria;
import fem2.mesh_and_topology.SurfaceSpatialCriteria;
import inf.math.LinAlg;
import inf.text.ArrayFormat;

import java.io.BufferedWriter;
import java.io.File;
import java.io.FileWriter;
import java.io.IOException;
import java.io.Writer;
import java.util.ArrayList;
import java.util.Collections;
import java.util.HashMap;
import java.util.HashSet;
import java.util.Map;
import java.util.Set;

import math2.FunctionRnToR;
import math2.MathUtilities;

/**
 * utility class to find nodes, edges, faces,... on a mesh
 * 
 * @author hbui
 * 
 */
public class MeshUtilities {

	/**
	 * find Node having coordinate (x, y, z)
	 * 
	 * @param x
	 * @return node
	 */
	public static Node seekNode(Mesh mesh, double... x) {
		for (int i = 0; i < mesh.countNodes(); i++) {
			Node n = mesh.getNode(i);
			double[] ncor = n.getPosition();
			double d = inf.math.LinAlg.diffNorm2(x.length, ncor, x);
			if (d < MathUtilities.tol14) {
				return n;
			}
		}
		return null;
	}

	/**
	 * find Node nearest to coordinate (x, y, z)
	 * 
	 * @param x
	 * @return node
	 */
	public static Node seekNodeNearest(Mesh mesh, double... x) {
		double d_old = Double.POSITIVE_INFINITY;
		Node n = null;
		for (int i = 0; i < mesh.countNodes(); i++) {
			Node n1 = mesh.getNode(i);
			double d = LinAlg.diffNorm2(x.length, n1.getPosition(), x);
			if (d < d_old) {
				n = n1;
				d_old = d;
			}
		}
		return n;
	}

	/**
	 * find node nearest to x but satisfied spatial criteria
	 * 
	 * @param mesh
	 * @param cr
	 * @param x
	 * @return
	 */
	public static Node seekNodeNearest(Mesh mesh, SpatialCriteria cr, double... x) {
		double d_old = Double.POSITIVE_INFINITY;
		Node n = null;
		for (int i = 0; i < mesh.countNodes(); i++) {
			Node n1 = mesh.getNode(i);
			double d = LinAlg.diffNorm2(x.length, n1.getPosition(), x);
			if (d < d_old && cr.isSatisfied(n1.getPosition())) {
				n = n1;
				d_old = d;
			}
		}
		return n;
	}

	/**
	 * 
	 * find Nodes on a surface describing by a function
	 * 
	 * @param f
	 * @return nodes
	 */
	public static ArrayList<Node> seekNodes(Mesh mesh, FunctionRnToR f) {
		ArrayList<Node> nodes = new ArrayList<Node>();
		for (int i = 0; i < mesh.countNodes(); i++) {
			Node n = mesh.getNode(i);
			double[] p = n.getPosition();
			double d = Math.abs(f.valueAt(p));
			if (d < MathUtilities.tol14) {
				nodes.add(n);
			}
		}
		return nodes;
	}

	/**
	 * 
	 * find Nodes on a surface describing by a spatial criteria
	 * 
	 * @param f
	 * @return nodes
	 */
	public static ArrayList<Node> seekNodes(Mesh mesh, SpatialCriteria c) {
		ArrayList<Node> nodes = new ArrayList<Node>();
		for (int i = 0; i < mesh.countNodes(); i++) {
			Node n = mesh.getNode(i);
			double[] p = n.getPosition();
			if (c.isSatisfied(p)) {
				nodes.add(n);
			}
		}
		return nodes;
	}

	/**
	 * find Node on surface (line in 2d) ax+by+c = 0
	 * 
	 * @param a
	 * @param b
	 * @param c
	 * @return list of Nodes
	 */
	public static ArrayList<Node> seekNodesOnSurface(Mesh mesh, double a, double b, double c) {
		// ArrayList<Node> nodes = new ArrayList<Node>();
		// for (int i = 0; i < mesh.countNodes(); i++) {
		// Node n = mesh.getNode(i);
		// double[] p = n.getPosition();
		// double f = Math.abs(a * p[0] + b * p[1] + c);
		// if (f < LinAlgUtilities.tol14) {
		// nodes.add(n);
		// }
		// }
		// return nodes;
		return seekNodes(mesh, new SurfaceSpatialCriteria(a, b, c));
	}

	/**
	 * find nodes on surface (line in 2d) ax+by+c = 0
	 * 
	 * @param a
	 * @param b
	 * @param c
	 * @return array of Nodes
	 */
	public static Node[] seekNodesOnSurfaceAsArray(Mesh mesh, double a, double b, double c) {
		ArrayList<Node> nodes = seekNodesOnSurface(mesh, a, b, c);
		Node[] nodesArray = new Node[nodes.size()];
		nodes.toArray(nodesArray);
		return nodesArray;
	}

	/**
	 * find nodes on surface ax+by+cz+d = 0
	 * 
	 * @param a
	 * @param b
	 * @param c
	 * @param d
	 * @return list of nodes
	 */
	public static ArrayList<Node> seekNodesOnSurface(Mesh mesh, double a, double b, double c,
			double d) {
		// ArrayList<Node> nodes = new ArrayList<Node>();
		// for (int i = 0; i < mesh.countNodes(); i++) {
		// Node n = mesh.getNode(i);
		// double[] p = n.getPosition();
		// double f = Math.abs(a * p[0] + b * p[1] + c * p[2] + d);
		// if (f < LinAlgUtilities.tol14) {
		// nodes.add(n);
		// }
		// }
		// return nodes;
		return seekNodes(mesh, new SurfaceSpatialCriteria(a, b, c, d));
	}

	/**
	 * find nodes on line cross (x0, y0) with direction (u, v)
	 * 
	 * @param mesh
	 * @param x0
	 * @param y0
	 * @param u
	 * @param v
	 * @return
	 */
	public static ArrayList<Node> seekNodesOnLine(Mesh mesh, double x0, double y0, double u,
			double v) {
		// ArrayList<Node> nodes = new ArrayList<Node>();
		// for (int i = 0; i < mesh.countNodes(); i++) {
		// Node n = mesh.getNode(i);
		// double[] p = n.getPosition();
		// double t;
		// if (Math.abs(u) > LinAlgUtilities.tol14) {
		// t = (p[0] - x0) / u;
		// } else {
		// if (Math.abs(v) > LinAlgUtilities.tol14) {
		// t = (p[1] - y0) / v;
		// } else {
		// throw new Error("direction is vector zero. How come?");
		// }
		// }
		// double dx = p[0] - (x0 + u * t);
		// double dy = p[1] - (y0 + v * t);
		// if (Math.sqrt(dx * dx + dy * dy) < LinAlgUtilities.tol14) {
		// nodes.add(n);
		// }
		// }
		// return nodes;
		return seekNodes(mesh, new LineSpatialCriteria(x0, y0, u, v));
	}

	/**
	 * find nodes on line cross (x0, y0, z0) with direction (u, v, w)
	 * 
	 * @param mesh
	 * @param x0
	 * @param y0
	 * @param z0
	 * @param u
	 * @param v
	 * @param w
	 * @return
	 */
	public static ArrayList<Node> seekNodesOnLine(Mesh mesh, double x0, double y0, double z0,
			double u, double v, double w) {
		// ArrayList<Node> nodes = new ArrayList<Node>();
		// for (int i = 0; i < mesh.countNodes(); i++) {
		// Node n = mesh.getNode(i);
		// double[] p = n.getPosition();
		// double t;
		// if (Math.abs(u) > LinAlgUtilities.tol14) {
		// t = (p[0] - x0) / u;
		// } else {
		// if (Math.abs(v) > LinAlgUtilities.tol14) {
		// t = (p[1] - y0) / v;
		// } else {
		// if (Math.abs(w) < LinAlgUtilities.tol14) {
		// throw new Error("direction is vector zero. How come?");
		// }
		// t = (p[2] - z0) / w;
		// }
		// }
		// double dx = p[0] - (x0 + u * t);
		// double dy = p[1] - (y0 + v * t);
		// double dz = p[2] - (z0 + w * t);
		// if (Math.sqrt(dx * dx + dy * dy + dz * dz) < LinAlgUtilities.tol14) {
		// nodes.add(n);
		// }
		// }
		// return nodes;
		return seekNodes(mesh, new LineSpatialCriteria(x0, y0, z0, u, v, w));
	}

	/**
	 * 
	 * @return list of independent nodes
	 */
	private static ArrayList<Integer> getIndependentNodes(Mesh mesh) {
		int[] nodesFlag = new int[mesh.countNodes()];
		/*
		 * check all blocks for nodes
		 */
		for (int i = 0; i < mesh.countBlocks(); i++) {
			Block b = mesh.getBlock(i);
			for (Node n : b.getNodes()) {
				nodesFlag[n.getId() - 1] = 1;
			}
		}
		/*
		 * check all faces for nodes
		 */
		for (int i = 0; i < mesh.countFaces(); i++) {
			Face f = mesh.getFace(i);
			for (Node n : f.getNodes()) {
				nodesFlag[n.getId() - 1] = 1;
			}
		}
		/*
		 * check all edges for nodes
		 */
		for (int i = 0; i < mesh.countEdges(); i++) {
			Edge e = mesh.getEdge(i);
			for (Node n : e.getNodes()) {
				nodesFlag[n.getId() - 1] = 1;
			}
		}
		/*
		 * extract independent nodes
		 */
		ArrayList<Integer> list = new ArrayList<Integer>();
		for (int i = 0; i < nodesFlag.length; i++) {
			if (nodesFlag[i] == 0) {
				list.add(i);
			}
		}
		return list;
	}

	/**
	 * check if mesh contains nodes not belong to and edges/faces/blocks
	 * 
	 * @return false if mesh is invalid<br>
	 *         true if valid
	 */
	public static boolean validateMesh(Mesh mesh) {
		ArrayList<Integer> list = getIndependentNodes(mesh);
		for (int i = 0; i < list.size(); i++) {
			Debugger.warn("Node " + (list.get(i) + 1) + " does not belong to any mesh part");
		}
		return list.size() > 0;
	}

	/**
	 * correct the mesh by removing independent nodes. Use it with care
	 */
	public static void correctMesh(Mesh mesh) {
		ArrayList<Integer> list = getIndependentNodes(mesh);
		Collections.sort(list);
		for (int i = list.size() - 1; i > -1; i--) {
			Debugger.warn("Remove node " + (list.get(i) + 1) + " from mesh");
			mesh.removeNode(list.get(i));
		}
	}

	/**
	 * mapping from MeshPartType to Gmsh element type<blockquote>refer to Gmsh
	 * manual, page 83</blockquote>
	 * 
	 * @param type
	 * @return elmType
	 */
	public static int getGmshElmType(MeshPartType type) {
		switch (type) {
		case LINE2:
			return 1;
		case LINE3:
			return 8;
		case TRI3:
			return 2;
		case TRI6:
			return 9;
		case QUAD4:
			return 3;
		case QUAD8:
			return 16;
		case TET4:
			return 4;
		case TET10:
			return 11;
		case HEX8:
			return 5;
		case HEX20:
			return 17;
		default:
			throw new Error("invalid mesh part type");
		}
	}

	/**
	 * mapping from MeshPartType to Comet element type
	 * 
	 * @param type
	 * @return string of element type
	 */
	public static String getCometElmType(MeshPartType type) {
		switch (type) {
		case LINE2:
			return "line";
		case LINE3:
			return "line";
		case TRI3:
			return "tri";
		case TRI6:
			return "tri";
		case QUAD4:
			return "quad";
		case QUAD8:
			return "quad";
		case TET4:
			return "tet";
		case TET10:
			return "tet";
		case HEX8:
			return "hex";
		case HEX20:
			return "hex";
		default:
			throw new Error("invalid mesh part type");
		}
	}

	/**
	 * write mesh to Gmsh format
	 * 
	 * @param mesh
	 * @param fileName
	 * @throws IOException
	 */
	public static void writeGmshMesh(Mesh mesh, String fileName) throws IOException {
		int nn = mesh.countNodes();
		int ne = mesh.countEdges();
		int nf = mesh.countFaces();
		int nb = mesh.countBlocks();

		File file = new File(fileName);
		int dim = mesh.getDimension();

		if (file.exists()) {
			/*
			 * delete mesh file
			 */
			file.delete();
		}
		/*
		 * create new file
		 */
		Writer output = new BufferedWriter(new FileWriter(file));
		/*
		 * write headers
		 */
		output.write("$MeshFormat\n");
		output.write("2.2 0 8\n");
		output.write("$EndMeshFormat\n");
		/*
		 * write nodes
		 */
		output.write("$Nodes\n");
		output.write(nn + "\n");
		for (int i = 0; i < nn; i++) {
			Node node = mesh.getNode(i);
			output.write(node.getId() + "");
			for (int j = 0; j < dim; j++) {
				double tmp = node.getPosition()[j];
				output.write("\t" + tmp);
			}
			if (dim == 2) {
				output.write("\t" + 0.0);
			}
			output.write("\n");
		}
		output.write("$EndNodes\n");

		/*
		 * write mesh parts
		 */
		output.write("$Elements\n");
		output.write((ne + nf + nb) + "\n");
		for (int i = 0; i < ne + nf + nb; i++) {
			MeshPart mp = null;

			if (i < ne) {
				mp = mesh.getEdge(i);
			} else if (i < ne + nf) {
				mp = mesh.getFace(i - ne);
			} else if (i < ne + nf + nb) {
				mp = mesh.getBlock(i - ne - nf);
			}

			int elmType = getGmshElmType(mp.getType());
			output.write((i + 1) + " ");
			output.write(elmType + " 0");
			Node[] nodes = mp.getNodes();
			for (int j = 0; j < nodes.length; j++) {
				output.write(" " + nodes[j].getId());
			}
			output.write("\n");
		}
		output.write("$EndElements\n");

		output.close();

		Debugger.watch("Write Gmsh mesh file complete.");
	}

	/**
	 * write a mesh to a format compatible with GiD. If you already have a
	 * model, use GidPost instead
	 * 
	 * @param mesh
	 * @param fileName
	 * @throws IOException
	 */
	public static void writeGiDMesh(Mesh mesh, String fileName) throws IOException {

		File file = new File(fileName);
		int dim = mesh.getDimension();

		if (file.exists()) {
			/*
			 * delete mesh file
			 */
			file.delete();
		}
		/*
		 * create new file
		 */
		Writer output = new BufferedWriter(new FileWriter(file));

		MeshPartType typicalType = null;
		Set<MeshPartType> typeSet = new HashSet<MeshPartType>();
		int ne;
		if (dim == 2) {
			ne = mesh.countFaces();
		} else if (dim == 3) {
			ne = mesh.countBlocks();
		} else {
			output.close();
			throw new Error("invalid dimension");
		}
		for (int i = 0; i < ne; i++) {
			MeshPart e = null;
			if (dim == 2) {
				e = mesh.getFace(i);
			} else if (dim == 3) {
				e = mesh.getBlock(i);
			}
			/*
			 * get the first mesh part type
			 */
			if (i == 0) {
				typicalType = e.getType();
			}
			/*
			 * check if there is other type of mesh part
			 */
			if (e.getType() != typicalType) {
				typeSet.add(e.getType());
			}

			if (i == 0) {
				/*
				 * write the header
				 */
				output.write("MESH dimension " + dim + " ElemType ");

				switch (typicalType) {
				case TRI3:
					output.write("Triangle Nnode 3\n");
					break;
				case TRI6:
					output.write("Triangle Nnode 6\n");
					break;
				case QUAD4:
					output.write("Quadrilateral Nnode 4\n");
					break;
				case QUAD8:
					output.write("Quadrilateral Nnode 8\n");
					break;
				case QUAD9:
					output.write("Quadrilateral Nnode 9\n");
					break;
				case TET4:
					output.write("Tetrahedra Nnode 4\n");
					break;
				case TET10:
					output.write("Tetrahedra Nnode 10\n");
					break;
				case HEX8:
					output.write("Hexahedra Nnode 8\n");
					break;
				case HEX20:
					output.write("Hexahedra Nnode 20\n");
					break;
				default:
					throw new Error("this type of mesh part is not supported. Is it new?");
				}

				/*
				 * write coordinates
				 */
				output.write("Coordinates\n");
				for (int j = 0; j < mesh.countNodes(); j++) {
					Node n = mesh.getNode(j);
					output.write(n.getId() + "");
					for (int k = 0; k < dim; k++) {
						output.write("\t" + n.getPosition()[k]);
					}
					output.write("\n");
				}
				output.write("End Coordinates\n");

				/*
				 * write typical elements
				 */
				output.write("Elements\n");
			}
			Node[] nodes = e.getNodes();
			if (e.getType() == typicalType) {
				output.write(e.getId() + "");
				for (int j = 0; j < e.countModes(); j++) {
					output.write("\t" + nodes[j].getId());
				}
				output.write("\n");
			}
		}
		output.write("End Elements\n");

		/*
		 * write other type of mesh part
		 */
		for (MeshPartType type : typeSet) {
			/*
			 * write the header
			 */
			output.write("MESH dimension " + dim + " ElemType ");

			switch (type) {
			case TRI3:
				output.write("Triangle Nnode 3\n");
				break;
			case TRI6:
				output.write("Triangle Nnode 6\n");
				break;
			case QUAD4:
				output.write("Quadrilateral Nnode 4\n");
				break;
			case QUAD8:
				output.write("Quadrilateral Nnode 8\n");
				break;
			case QUAD9:
				output.write("Quadrilateral Nnode 9\n");
				break;
			case TET4:
				output.write("Tetrahedra Nnode 4\n");
				break;
			case TET10:
				output.write("Tetrahedra Nnode 10\n");
				break;
			case HEX8:
				output.write("Hexahedra Nnode 8\n");
				break;
			case HEX20:
				output.write("Hexahedra Nnode 20\n");
				break;
			default:
				throw new Error("this type of mesh part is not supported. Is it new?");
			}

			for (int i = 0; i < ne; i++) {
				MeshPart e = null;
				if (dim == 2) {
					e = mesh.getFace(i);
				} else if (dim == 3) {
					e = mesh.getBlock(i);
				}
				/*
				 * check type of mesh part
				 */
				if (e.getType() == type) {
					/*
					 * write elements
					 */
					Node[] nodes = e.getNodes();
					output.write("Elements\n");
					if (e.getType() == typicalType) {
						output.write(e.getId() + "");
						for (int j = 0; j < e.countModes(); j++) {
							output.write("\t" + nodes[j].getId());
						}
						output.write("\n");
					}
				}

			}
			output.write("End Elements\n");
		}

		/*
		 * close
		 */
		output.close();

		Debugger.watch("Write GiD mesh file complete.");
	}

	/**
	 * write mesh to Comet format
	 * 
	 * @param mesh
	 * @param fileName
	 * @throws IOException
	 */
	public static void writeCometMesh(Mesh mesh, String fileName) throws IOException {
		int nn = mesh.countNodes();
		int nf = mesh.countFaces();
		int nb = mesh.countBlocks();

		File file = new File(fileName);
		int dim = mesh.getDimension();

		if (file.exists()) {
			/*
			 * delete mesh file
			 */
			file.delete();
		}
		/*
		 * create new file
		 */
		Writer output = new BufferedWriter(new FileWriter(file));
		/*
		 * write headers
		 */
		if (dim == 2) {
			output.write(nn + " " + nf);
		} else if (dim == 3) {
			output.write(nn + " " + nb);
		}
		output.write(" 0 0 0\n");
		/*
		 * write nodes
		 */
		for (int i = 0; i < nn; i++) {
			Node node = mesh.getNode(i);
			output.write(node.getId() + "");
			for (int j = 0; j < dim; j++) {
				double tmp = node.getPosition()[j];
				output.write("\t" + tmp);
			}
			if (dim == 2) {
				output.write("\t" + 0.0);
			}
			output.write("\n");
		}
		/*
		 * write mesh parts
		 */
		int ne;
		if (dim == 2) {
			ne = nf;
		} else if (dim == 3) {
			ne = nb;
		} else {
			output.close();
			throw new Error("invalid dimension");
		}
		for (int i = 0; i < ne; i++) {
			MeshPart mp = null;

			if (dim == 2) {
				mp = mesh.getFace(i);
			} else if (dim == 3) {
				mp = mesh.getBlock(i);
			}

			String elmType = getCometElmType(mp.getType());
			output.write((i + 1) + " 1 ");
			output.write(elmType);
			Node[] nodes = mp.getNodes();
			for (int j = 0; j < nodes.length; j++) {
				output.write(" " + nodes[j].getId());
			}
			output.write("\n");
		}

		output.close();

		Debugger.watch("Write Comet mesh file complete.");
	}

	/**
	 * join the 2 meshes at the boundary defined by a spatial criteria
	 * 
	 * @param mesh1
	 * @param mesh2
	 * @return
	 */
	public static Mesh joinMesh(Mesh mesh1, Mesh mesh2, SpatialCriteria c) {
		/*
		 * find nodes on mesh1 on boundary f
		 */
		ArrayList<Node> b1 = seekNodes(mesh1, c);
		/*
		 * find nodes on mesh2 on boundary f
		 */
		ArrayList<Node> b2 = seekNodes(mesh2, c);
		/*
		 * create a map from b2 to b1
		 */
		Map<Node, Node> mapNodes = new HashMap<Node, Node>();
		for (int i = 0; i < b2.size(); i++) {
			Node n = b2.get(i);
			int j;
			for (j = 0; j < b1.size(); j++) {
				Node n1 = b1.get(j);
				if (n.equals(n1)) {
					mapNodes.put(n, n1);
					break;
				}
			}
			if (j == b1.size()) {
				throw new Error("unmatched node on mesh2 found at"
						+ ArrayFormat.format(n.getPosition()));
			}
		}
		/*
		 * set the mesh to mesh1
		 */
		Mesh mesh = mesh1;
		/*
		 * add all nodes from mesh2 which does not belong to b2
		 */
		for (int i = 0; i < mesh2.countNodes(); i++) {
			Node n = mesh2.getNode(i);
			if (!(b2.contains(n))) {
				mesh.addNode(n);
			}
		}
		/*
		 * get all edges from mesh2, replace nodes and add to mesh1
		 */
		for (int i = 0; i < mesh2.countEdges(); i++) {
			Edge e = mesh2.getEdge(i);
			Node[] nodes = e.getNodes();
			for (int j = 0; j < nodes.length; j++) {
				if (b2.contains(nodes[j])) {
					/*
					 * replace by according node in mesh1
					 */
					nodes[j] = mapNodes.get(nodes[j]);
				}
			}
			mesh.addEdge(e);
		}
		/*
		 * get all faces from mesh2, replace nodes and add to mesh1
		 */
		for (int i = 0; i < mesh2.countFaces(); i++) {
			Face f = mesh2.getFace(i);
			Node[] nodes = f.getNodes();
			for (int j = 0; j < nodes.length; j++) {
				if (b2.contains(nodes[j])) {
					/*
					 * replace by according node in mesh1
					 */
					nodes[j] = mapNodes.get(nodes[j]);
				}
			}
			mesh.addFace(f);
		}
		/*
		 * get all blocks from mesh2, replace nodes and add to mesh1
		 */
		for (int i = 0; i < mesh2.countBlocks(); i++) {
			Block b = mesh2.getBlock(i);
			Node[] nodes = b.getNodes();
			for (int j = 0; j < nodes.length; j++) {
				if (b2.contains(nodes[j])) {
					/*
					 * replace by according node in mesh1
					 */
					nodes[j] = mapNodes.get(nodes[j]);
				}
			}
			mesh.addBlock(b);
		}

		return mesh;
	}

	/**
	 * create list of edges from list of nodes. Nodes need to be sorted before
	 * creating edges. The number of fields of nodes must be consistent.
	 * 
	 * @param nodes
	 * @return
	 */
	public static ArrayList<Edge> createEdges(ArrayList<Node> nodes) {
		ArrayList<Edge> edges = new ArrayList<Edge>();
		int d = nodes.get(0).countFields();

		for (int i = 0; i < nodes.size() - 1; i++) {
			Node n1 = nodes.get(i);
			Node n2 = nodes.get(i + 1);
			Edge e = new Edge(n1, n2);
			e.setNumberOfFields(d);
			edges.add(e);
		}

		return edges;
	}

}
